! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

MODULE CONV_FORCING_MODS

CONTAINS

SUBROUTINE GET_POINTS_INFO(INFO_FILE,CHECK,NPOINTS,PPOINTS,LAT,LON)
USE NETCDF
IMPLICIT NONE
CHARACTER(LEN=*), INTENT(IN) :: INFO_FILE
LOGICAL,INTENT(IN)           :: CHECK
INTEGER,INTENT(OUT)          :: NPOINTS
INTEGER,DIMENSION(:),ALLOCATABLE,INTENT(OUT)      :: PPOINTS
REAL,DIMENSION(:),ALLOCATABLE,INTENT(OUT) :: LAT,LON


INTEGER :: NCID,VARID,DIMID
INTEGER :: STATUS

REAL,DIMENSION(:),ALLOCATABLE :: LSM,CLAKE
INTEGER                               :: JP,NPLAND,NPLAKE,NPLAKE_FRAC,NPOCEAN
INTEGER                               :: LATMIN,LONMIN,LATMAX,LONMAX

!*==============================================0
!* OPEN INFO FILE
CALL NCERROR( NF90_OPEN(INFO_FILE,NF90_NOWRITE,NCID), 'OPENING INFO FILE SUB_GPI' )
!* GET X DIMENSION LENGTH - NPOINTS ! 
CALL NCERROR( NF90_INQ_DIMID(NCID,'x',DIMID), 'GETTING x DIMID INFO FILE ')
CALL NCERROR( NF90_INQUIRE_DIMENSION(ncid=NCID,dimid=DIMID,len=NPOINTS), 'GETTING x LEN INFO FILE')

!* ALLOCATE STUFF
IF( .NOT. ALLOCATED(PPOINTS)) ALLOCATE(PPOINTS(NPOINTS))
IF( .NOT. ALLOCATED(LAT))     ALLOCATE(LAT(NPOINTS))
IF( .NOT. ALLOCATED(LON))     ALLOCATE(LON(NPOINTS))

!* GET PPOINTS, LAT AND LON_TOT
CALL NCERROR(NF90_INQ_VARID(NCID,'x',VARID),'GETTING x VARID INFO FILE')
CALL NCERROR(NF90_GET_VAR(NCID,VARID,PPOINTS(1:NPOINTS)), 'READING x VAR INFO FILE')
CALL NCERROR(NF90_INQ_VARID(NCID,'lat',VARID),'GETTING lat VARID INFO FILE')
CALL NCERROR(NF90_GET_VAR(NCID,VARID,LAT(1:NPOINTS)), 'READING lat VAR INFO FILE')
CALL NCERROR(NF90_INQ_VARID(NCID,'lon',VARID),'GETTING lon VARID INFO FILE')
CALL NCERROR(NF90_GET_VAR(NCID,VARID,LON(1:NPOINTS)), 'READING LON VAR INFO FILE')


IF ( CHECK ) THEN
  PRINT*, "=====CHECK IN GET_POINTS_INFO====="
  PRINT*,'NPOINTS=',NPOINTS
  !* GET LSM AND CLAKE
  ALLOCATE(LSM(NPOINTS))
  ALLOCATE(CLAKE(NPOINTS))
  CALL NCERROR(NF90_INQ_VARID(NCID,'landsea',VARID),'GETTING landsea VARID INFO FILE')
  CALL NCERROR(NF90_GET_VAR(NCID,VARID,LSM(1:NPOINTS)), 'READING lsm VAR INFO FILE')
  CALL NCERROR(NF90_INQ_VARID(NCID,'CLAKE',VARID),'GETTING CLAKE VARID INFO FILE')
  CALL NCERROR(NF90_GET_VAR(NCID,VARID,CLAKE(1:NPOINTS)), 'READING CLAKE VAR INFO FILE')
  !* CALC NP*
  NPLAND=0
  NPLAKE=0
  NPLAKE_FRAC=0
  NPOCEAN=0
  LATMAX=LAT(1)
  LONMAX=LON(1)
  LATMIN=LAT(1)
  LONMIN=LON(1)
  DO JP=1,NPOINTS
    !PRINT*,LSM(JP),CLAKE(JP),PPOINTS(JP),LAT(JP),LON(JP)
    IF (LSM(JP) > 0.5) THEN
      NPLAND=NPLAND+1
      IF( CLAKE(JP) > 0.5*(1.-LSM(JP)) ) THEN
        NPLAKE_FRAC=NPLAKE_FRAC+1
      ENDIF
    ELSEIF( CLAKE(JP) > 0.5 ) THEN
      NPLAKE=NPLAKE+1
    ELSE
      NPOCEAN=NPOCEAN+1
    ENDIF
    IF ( LAT(JP) < LATMIN ) LATMIN=LAT(JP)
    IF ( LAT(JP) > LATMAX ) LATMAX=LAT(JP)
    IF ( LONMIN < LONMAX ) THEN
      IF ( LON(JP) < LONMIN ) LONMIN=LON(JP)
      IF ( LON(JP) > LONMAX ) LONMAX=LON(JP)
    ELSE
      IF ( LON(JP) < LONMIN .AND. LON(JP).GT.180. ) LONMIN=LON(JP)
      IF ( LON(JP) > LONMAX .AND. LON(JP).LT.180. ) LONMAX=LON(JP)
    ENDIF
  ENDDO
  
  PRINT*,'NPLAND=',NPLAND
  PRINT*,'NPLAKE=',NPLAKE
  PRINT*,'NPLAKE_FRAC=',NPLAKE_FRAC
  PRINT*,'NPOCEAN',NPOCEAN
  PRINT*,'LATMIN=',LATMIN
  PRINT*,'LATMAX=',LATMAX
  PRINT*,'LONMIN=',LONMIN
  PRINT*,'LONMAX=',LONMAX
  PRINT*, "=====CHECK IN GET_POINTS_INFO====="
  DEALLOCATE(LSM)
  DEALLOCATE(CLAKE)
ENDIF

CALL NCERROR( NF90_CLOSE(NCID), 'CLOSING INFO FILE SUB_GPI')

END SUBROUTINE GET_POINTS_INFO


SUBROUTINE NCERROR(STATUS,STRING)
USE NETCDF
IMPLICIT NONE

INTEGER,INTENT(IN) :: STATUS
CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: STRING 

IF ( STATUS /= 0 ) THEN 
  PRINT*, TRIM(NF90_STRERROR(STATUS))
  IF( PRESENT(STRING) ) PRINT*,TRIM(STRING)
  PRINT*,'PROGRAM STOP ! '
  STOP 9
ENDIF

END SUBROUTINE NCERROR

END MODULE CONV_FORCING_MODS



PROGRAM CONV_FORCING
USE NETCDF
USE GRIB_API
USE CONV_FORCING_MODS
IMPLICIT NONE 

CHARACTER(LEN=500)                       :: GRIB_FILE,NETCDF_FILE,VAR_NAME,INFO_FILE
INTEGER                                  :: NPOINTS
INTEGER,DIMENSION(:),ALLOCATABLE         :: PPOINTS
REAL,DIMENSION(:),ALLOCATABLE            :: LAT,LON

INTEGER                                  :: STATUS,NCID,VARID,DIMID
INTEGER                                  :: IGRIB,ERR,IFILE,NPROD
INTEGER                                  :: TSTEP,NTSTEPS
REAL                                     :: TIME,ZDTFORC
REAL,DIMENSION(:),ALLOCATABLE            :: VALUES,DATA_VALUES
INTEGER                                  :: YEAR,MONTH,DAY,HOUR,STEP,NPOINTS_TOT
INTEGER                                  :: YYYYMM,MONTHS1(12),MONTHS2(12),YYEAR,MMONTH

INTEGER :: JPP,JP

DATA MONTHS2/31,29,31,30,31,30,31,31,30,31,30,31/
DATA MONTHS1/31,28,31,30,31,30,31,31,30,31,30,31/


NAMELIST /INPUT/GRIB_FILE,NETCDF_FILE,VAR_NAME,INFO_FILE,ZDTFORC


!* ======================================================
!* OPEN NAMELIST FILE 
OPEN(10,FILE='input.nam',STATUS='OLD')
READ(10,INPUT)
CLOSE(10)
 
PRINT*,'READING INFO FROM NETCDF FILE'
PRINT*,TRIM(INFO_FILE)
PRINT*,'READING FROM GRIB FILE'
PRINT*, TRIM(GRIB_FILE)
PRINT*,'WRITING TO NETCDF FILE'
PRINT*, TRIM(NETCDF_FILE)


!* GET INFO ON POINTS TO GET FROM GRIB FILE 
!=============================================

CALL GET_POINTS_INFO(INFO_FILE,.TRUE.,NPOINTS,PPOINTS,LAT,LON)
 
 
!* ===========================================================
!* OPEN NETCDF FILE 
CALL NCERROR( NF90_OPEN(TRIM(NETCDF_FILE),NF90_WRITE,NCID),'OPEN NETCDF OUTPUT' )
CALL NCERROR( NF90_INQ_DIMID(NCID,'x',DIMID),'GET land DIMID' )
CALL NCERROR( NF90_INQUIRE_DIMENSION(NCID,DIMID,LEN=JPP),'GET land LEN' )
IF ( JPP /= NPOINTS ) THEN
PRINT*, 'NUMBER OF POINTS IN NETCDF=',JPP, 'AND FROM GRIB=',NPOINTS
PRINT*, 'THEY MUST BE EQUAL!!! EXITING'
STOP 9 
ENDIF
! 
!* START WRITING LAT,LON,NPOINTS TO THE FILE
CALL NCERROR( NF90_INQ_VARID(NCID,'lat',VARID),'GET lat ID' )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,LAT,(/1/),(/NPOINTS/)),'WRITE lat' )
CALL NCERROR( NF90_INQ_VARID(NCID,'lon',VARID),'GET lon ID' )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,LON,(/1/),(/NPOINTS/)),'WRITE lon' )
CALL NCERROR( NF90_INQ_VARID(NCID,'x',VARID),'GET land ID' )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,PPOINTS,(/1/),(/NPOINTS/)),'WRITE x' )
! 
! 
!* ===============================================================
!* OPEN MAIN GRIB FILE AND DO THINGS ....

ALLOCATE(DATA_VALUES(NPOINTS))
IFILE=1
CALL GRIB_OPEN_FILE(IFILE,TRIM(GRIB_FILE),'R')
NPROD=0
TSTEP=0
CALL GRIB_NEW_FROM_FILE(IFILE,IGRIB,ERR)
DO WHILE (ERR /= GRIB_END_OF_FILE )
  NPROD=NPROD+1
  CALL GRIB_GET(IGRIB,'year',YEAR)
  CALL GRIB_GET(IGRIB,'month',MONTH)
  CALL GRIB_GET(IGRIB,'day',DAY)
  CALL GRIB_GET(IGRIB,'hour',HOUR)
  CALL GRIB_GET(IGRIB,'step',STEP)
  CALL GRIB_GET(IGRIB,'numberOfValues',NPOINTS_TOT)
  IF ( NPOINTS_TOT < NPOINTS ) THEN
    PRINT*,'NPOINTS_TOT < NPOINTS ->:',NPOINTS_TOT,NPOINTS
    PRINT*,'STOPING! '
    STOP 10
  ENDIF
  
  IF ( .NOT. ALLOCATED(VALUES) ) ALLOCATE(VALUES(NPOINTS_TOT))
  CALL GRIB_GET(IGRIB,'values',VALUES)

  SELECT CASE (TRIM(VAR_NAME))
    CASE ('SWdown','LWdown','Snowf','Rainf','Tair','CO2air','Qair','Wind_E','Wind_N', &
    & 'Wind','Ctpf','Ascat00', 'Ascat06', 'Ascat12', 'Ascat18', 'Ascat001', &
    & 'H_Ascat00','H_Ascat06', 'H_Ascat12','H_Ascat18','H_Ascat001','AWC', &
    & 'T2m00', 'T2m06', 'T2m12', 'T2m18', 'RH2m00', 'RH2m06', 'RH2m12', 'RH2m18', &
    & 'N_Ascat00', 'N_Ascat06', 'N_Ascat12', 'N_Ascat18', 'N_Ascat001', &
    & 'sekf_ascat_bias_in', &
    & 'covar_ssm_00_swvl1', 'covar_ssm_00_swvl2', 'covar_ssm_00_swvl3', &
    & 'covar_ssm_06_swvl1', 'covar_ssm_06_swvl2', 'covar_ssm_06_swvl3', &
    & 'covar_ssm_12_swvl1', 'covar_ssm_12_swvl2', 'covar_ssm_12_swvl3', &
    & 'covar_ssm_18_swvl1', 'covar_ssm_18_swvl2', 'covar_ssm_18_swvl3', &
    & 'var_00_swvl1', 'var_00_swvl2', 'var_00_swvl3', &
    & 'var_06_swvl1', 'var_06_swvl2', 'var_06_swvl3', &
    & 'var_12_swvl1', 'var_12_swvl2', 'var_12_swvl3', &
    & 'var_18_swvl1', 'var_18_swvl2', 'var_18_swvl3', &
    & 'covar_T2m_00_swvl1', 'covar_T2m_00_swvl2', 'covar_T2m_00_swvl3', &
    & 'covar_T2m_06_swvl1', 'covar_T2m_06_swvl2', 'covar_T2m_06_swvl3', &
    & 'covar_T2m_12_swvl1', 'covar_T2m_12_swvl2', 'covar_T2m_12_swvl3', &
    & 'covar_T2m_18_swvl1', 'covar_T2m_18_swvl2', 'covar_T2m_18_swvl3', &
    & 'covar_RH2m_00_swvl1', 'covar_RH2m_00_swvl2', 'covar_RH2m_00_swvl3', &
    & 'covar_RH2m_06_swvl1', 'covar_RH2m_06_swvl2', 'covar_RH2m_06_swvl3', &
    & 'covar_RH2m_12_swvl1', 'covar_RH2m_12_swvl2', 'covar_RH2m_12_swvl3', &
    & 'covar_RH2m_18_swvl1', 'covar_RH2m_18_swvl2', 'covar_RH2m_18_swvl3')
      TSTEP=TSTEP+1
    CASE ('PSurf')
      TSTEP=TSTEP+1
      VALUES(1:NPOINTS_TOT)=EXP(VALUES(1:NPOINTS_TOT))
      IF ( TSTEP .EQ. 1 ) PRINT*,'--> PSurf=EXP(PSurf) <--'
    CASE DEFAULT
      PRINT*,'VARIABLE NOT LISTEDD!!!!'
      STOP 9
  END SELECT 
  TIME=(TSTEP-1)*ZDTFORC 
  DO JP=1,NPOINTS
    DATA_VALUES(JP)=VALUES(PPOINTS(JP))
  ENDDO

  CALL NCERROR( NF90_INQ_VARID(NCID,'time',VARID), 'get timeId')
  CALL NCERROR( NF90_PUT_VAR(NCID,VARID,TIME,(/TSTEP/)),'write time' ) 

  CALL NCERROR( NF90_INQ_VARID(NCID,TRIM(VAR_NAME),VARID),'get vId' )
  CALL NCERROR( NF90_PUT_VAR(NCID,VARID,DATA_VALUES,(/1,TSTEP/),(/NPOINTS,1/)),'write varId' )

  CALL GRIB_RELEASE(IGRIB)
  CALL GRIB_NEW_FROM_FILE(IFILE,IGRIB,ERR) 
  WRITE(*,'(7(I6,1X),F12.1)')YEAR,MONTH,DAY,HOUR,STEP,NPOINTS_TOT,TSTEP,TIME  
ENDDO 

PRINT*,NPOINTS,'NPOINTS CONVERTED OVER ',TSTEP,'TIME STEPS'
PRINT*,NPROD,'GRIB FIELDS DECODED'

CALL GRIB_CLOSE_FILE(IFILE)
CALL NCERROR( NF90_CLOSE(NCID),'CLOSE NETCDF FILE' )
 


END PROGRAM CONV_FORCING 
